Goal

Compare eQTLs with BXD-specific references and B6mm10 reference:

Information

Doing this in default STAR setting (genotypesandimputed_withannotation_Local_0_10) because also available with B6 mm10 reference (abbreviated B6)

FastQTL output (http://fastqtl.sourceforge.net/):

Preparation

Functions definitions

prepareRawCounts <- function(counts, ref){
  
  # convert gene ids to gene names
  idx_uniq <- which(!duplicated(GeneConvert$V2[match(rownames(counts), GeneConvert$V1)]))
  counts <- counts[idx_uniq,]
  rownames(counts) <- GeneConvert$V2[match(rownames(counts), GeneConvert$V1)]
  
  # change X into C for cortex
  colnames(counts) <- gsub("X", "C", colnames(counts))
  # split by tissue and condition
  CortexNSD <- subset(counts, select=grep("^C[0-9]*nsd$", colnames(counts), value=TRUE))
  CortexSD <- subset(counts, select=grep("^C[0-9]*sd$", colnames(counts), value=TRUE))
  LiverNSD <- subset(counts, select=grep("^L[0-9]*nsd$", colnames(counts), value=TRUE))
  LiverSD <- subset(counts, select=grep("^L[0-9]*sd$", colnames(counts), value=TRUE))
  # clean and formate names of mouse lines
  colnames(CortexNSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(CortexNSD)), LinesNames$V10)]
  colnames(CortexSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(CortexSD)), LinesNames$V10)]
  colnames(LiverNSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(LiverNSD)), LinesNames$V10)]
  colnames(LiverSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(LiverSD)), LinesNames$V10)]
  
  # output 4 datasets
  assign(paste0("counts", ref, "ref_CortexNSD"), CortexNSD, envir=parent.frame())
  assign(paste0("counts", ref, "ref_CortexSD"), CortexSD, envir=parent.frame())
  assign(paste0("counts", ref, "ref_LiverNSD"), LiverNSD, envir=parent.frame())
  assign(paste0("counts", ref, "ref_LiverSD"), LiverSD, envir=parent.frame())
}

# find top 3 genes for affected
top3genes <- function(tissue, condition){
  t <- unlist(strsplit(tissue, split=""))[1]
  
  dataB6 <- eval(parse(text=(paste0("dataB6", t, condition))))
  dataB6full <- eval(parse(text=(paste0("dataB6full", t, condition))))
  dataBXD <- eval(parse(text=(paste0("dataBXD", t, condition))))
  dataBXDfull <- eval(parse(text=(paste0("dataBXDfull", t, condition))))

  # identify genes in common between the B6 and BXD lists
  commoneGenes <- intersect(rownames(dataB6), rownames(dataBXD))

  # remove genes with 0 variants
  commoneGenes <- setdiff(commoneGenes, rownames(dataB6full)[dataB6full$V2==0])

  # bigger difference in slope
  orderedgenes <- commoneGenes[order(abs(dataB6full[commoneGenes, "V9"]-dataBXDfull[commoneGenes, "V9"]), decreasing=TRUE)]
  ##dataB6full[match(orderedgenes[1:3], rownames(dataB6full)),]
  ##dataBXDfull[match(orderedgenes[1:3], rownames(dataBXDfull)),]
  outtable <- cbind(dataB6full[orderedgenes[1:3], c(5, 8)], dataB6[orderedgenes[1:3], "adjustedpvalue"], dataBXDfull[orderedgenes[1:3], c(5, 8)], dataBXD[orderedgenes[1:3], "adjustedpvalue"])
  colnames(outtable) <- c("variantB6", "slopeB6", "qvalueB6", "variantBXD", "slopeBXD", "qvalueBXD")
  
  # output a table with info on 3 top affected genes
  assign(paste0("top3genes", t, condition), outtable, envir=parent.frame())
}

# plot gene expression for a pair gene-variant
plotGeneExpressionVariant <- function(ref, gene, variant, tissue, condition){
  # get first letter of tissue
   t <- unlist(strsplit(tissue, split=""))[1]
  # get gene expression dataset
  expr <- eval(parse(text=paste0("expr", ref, "ref", t, condition)))
  # plot
  plot(as.numeric(expr[gene, sort(names(expr))]), pch=19, col=as.factor(genotypesGN[genotypesGN$Locus==variant, sort(names(expr))]), main=paste0("Gene: ", gene,", ref: ", ref, ", variant: ", variant, ", tissue: ", tissue, ", condition: ", condition), xaxt="n", xlab="", ylab="gene expression", las=1)
  axis(1, 1:length(names(expr)), sort(names(expr)), las=2)
}
# vectorize function
plotGeneExpressionVariant <- Vectorize(plotGeneExpressionVariant, SIMPLIFY=FALSE)

# plot raw gene counts for a pair gene-variant
plotGeneCountsVariant <- function(ref, gene, variant, tissue, condition){
  # get gene expression dataset
  counts <- eval(parse(text=paste0("counts", ref, "ref_", tissue, condition)))
  # plot
  plot(as.numeric(counts[gene, sort(names(counts))]), pch=19, col=as.factor(genotypesGN[genotypesGN$Locus==variant, sort(names(counts))]), main=paste0("Gene: ", gene,", ref: ", ref, ", variant: ", variant, ", tissue: ", tissue, ", condition: ", condition), xaxt="n", xlab="", ylab="raw gene counts", las=1)
  axis(1, 1:length(names(counts)), sort(names(counts)), las=2)
}
# vectorize function
plotGeneCountsVariant <- Vectorize(plotGeneCountsVariant, SIMPLIFY=FALSE)

Data preparation

# load table to convert mouse line names
LinesNames <- read.table("F:/BXD/data/ConvertLineNames.tsv", header=FALSE, stringsAsFactors=FALSE, sep="\t")

# load table to convert between gene id and gene name
GeneConvert <- read.table("F:/BXD/references/gene_names_B6mm10.tab", header=FALSE, stringsAsFactors=FALSE, sep="\t")

# retrieve GeneNetwork (GN) genotypes
genotypesGN <- read.table("F:/BXD/data/genome/BXDGenotypes.geno", header=TRUE)

# retrieve eQTL results
for(tissue in c("Cortex", "Liver")){
  t <- unlist(strsplit(tissue, split=""))[1]
  for(condition in c("NSD", "SD")){
    assign(paste0("dataB6", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition, "_pvalcorrected.txt"), row.names=1, stringsAsFactors=FALSE))
    assign(paste0("dataB6full", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition, ".txt"), row.names=1, stringsAsFactors=FALSE))
    assign(paste0("dataBXD", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition,"_pvalcorrected.txt"), row.names=1, stringsAsFactors=FALSE))
    assign(paste0("dataBXDfull", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition, ".txt"), row.names=1, stringsAsFactors=FALSE))
  }
}

# retrieve raw gene counts, formate and subset
countsB6ref <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/Summary_ReadsPerGene.out.tab", header=TRUE, row.names=1)
prepareRawCounts(countsB6ref, "B6")
countsBXDref <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/Summary_ReadsPerGene.out.tab", header=TRUE, row.names=1)
prepareRawCounts(countsBXDref, "BXD")

# retrieve gene expression
exprB6refCNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Cortex_NSD.tab")
exprBXDrefCNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Cortex_NSD.tab")
exprB6refCSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Cortex_SD.tab")
exprBXDrefCSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Cortex_SD.tab")
exprB6refLNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Liver_NSD.tab")
exprBXDrefLNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Liver_NSD.tab")
exprB6refLSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Liver_SD.tab")
exprBXDrefLSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Liver_SD.tab")

Quantifying

Cortex NSD

## [1] "Same number of genes"
## [1] 16382

## [1] 4624
## [1] 4631
## [1] 4627.5
## [1] 4431
## [1] 95.75365
## [1] 15807
## [1] 16382
## [1] 96.49005
## [1] 4334
## [1] 93.65748
## [1] 193
## [1] 200

## [1] 4323
## [1] 4431
## [1] 97.56263
## [1] 15745
## [1] 16382
## [1] 96.11159

## [1] 28.07294
## [1] 28.22664
## 
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
## 
## pi0: 0.5596393   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     2376   2988  3950   4629  5298 6321 16452
## q-value     2067   2614  3466   4000  4634 5539 16452
## local FDR   1668   2041  2613   2935  3231 3635 16440

## 
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
## 
## pi0: 0.5580651   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     2353   2957  3950   4629  5276 6282 16395
## q-value     2061   2594  3463   3999  4643 5501 16395
## local FDR   1635   2031  2599   2922  3250 3640 16379

Cortex SD

## [1] "Same number of genes"
## [1] 16382

## [1] 5018
## [1] 5020
## [1] 5019
## [1] 4834
## [1] 96.31401
## [1] 15729
## [1] 16382
## [1] 96.01392
## [1] 4727
## [1] 94.18211
## [1] 184
## [1] 186

## [1] 4714
## [1] 4834
## [1] 97.51758
## [1] 15667
## [1] 16382
## [1] 95.63545

## [1] 30.4598
## [1] 30.58545
## 
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
## 
## pi0: 0.5149926   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     2569   3195  4224   4882  5627 6577 16451
## q-value     2278   2862  3775   4398  5028 6054 16452
## local FDR   1812   2243  2853   3204  3518 3985 16443

## 
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
## 
## pi0: 0.5160343   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     2540   3196  4218   4887  5545 6555 16395
## q-value     2265   2852  3766   4411  5031 6012 16395
## local FDR   1798   2233  2851   3210  3526 3973 16388

Liver NSD

## [1] "Same number of genes"
## [1] 12752

## [1] 3192
## [1] 3202
## [1] 3197
## [1] 3077
## [1] 96.24648
## [1] 12467
## [1] 12752
## [1] 97.76506
## [1] 3032
## [1] 94.83891
## [1] 115
## [1] 125

## [1] 3026
## [1] 3077
## [1] 98.34254
## [1] 12435
## [1] 12752
## [1] 97.51412

## [1] 24.97083
## [1] 25.09767
## 
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
## 
## pi0: 0.5610873   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     1601   2031  2803   3289  3842 4647 12819
## q-value     1368   1743  2347   2783  3210 3905 12819
## local FDR   1102   1356  1751   1983  2207 2531 12795

## 
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
## 
## pi0: 0.5532508   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     1589   2019  2782   3283  3833 4642 12764
## q-value     1351   1735  2345   2775  3212 3941 12764
## local FDR   1074   1339  1744   1977  2202 2538 12739

Liver SD

## [1] "Same number of genes"
## [1] 12752

## [1] 2817
## [1] 2838
## [1] 2827.5
## [1] 2712
## [1] 95.91512
## [1] 12470
## [1] 12752
## [1] 97.78858
## [1] 2681
## [1] 94.81874
## [1] 105
## [1] 126

## [1] 2678
## [1] 2712
## [1] 98.74631
## [1] 12444
## [1] 12752
## [1] 97.58469

## [1] 22.0459
## [1] 22.24566
## 
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
## 
## pi0: 0.6395766   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     1361   1790  2540   3049  3566 4348 12819
## q-value     1119   1460  2015   2392  2834 3407 12819
## local FDR    905   1109  1475   1704  1903 2191 12794

## 
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
## 
## pi0: 0.6250782   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1    <1
## p-value     1364   1785  2538   3053  3547 4349 12764
## q-value     1119   1444  2032   2398  2847 3418 12764
## local FDR    893   1110  1467   1709  1919 2207 12747

compare with Data descriptor

## [1] "ceCNSD_mm9"
## [1] 14889
## [1] 4522
## [1] 30.37142
## [1] "ceCSD_mm9"
## [1] 14889
## [1] 4732
## [1] 31.78185
## [1] "ceLNSD_mm9"
## [1] 14103
## [1] 3155
## [1] 22.37113
## [1] "ceLSD_mm9"
## [1] 14103
## [1] 2654
## [1] 18.81869
## [1] "ceCNSD_mm10"
## [1] 15734
## [1] 4192
## [1] 26.64294
## [1] "ceCSD_mm10"
## [1] 15734
## [1] 4542
## [1] 28.86742
## [1] "ceLNSD_mm10"
## [1] 12647
## [1] 3092
## [1] 24.44849
## [1] "ceLSD_mm10"
## [1] 12647
## [1] 2695
## [1] 21.3094

Examples of genes

Legend

  • black = B6 allele
  • tan = D2 allele
  • red = heterozygous

Cortex NSD

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

Cortex SD

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

About 2% of the genes are highly affected. Knowing we affected less than 1% of bp in the genome.

legend: black=B6 allele, tan=D2 allele

Checking for raw gene counts. There are comparable since we used the same reads on B6 or BXD references.

Liver NSD

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

Liver SD

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $B6
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

## $BXD
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## 
## $<NA>
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33

Session information

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=French_Switzerland.1252  LC_CTYPE=French_Switzerland.1252   
## [3] LC_MONETARY=French_Switzerland.1252 LC_NUMERIC=C                       
## [5] LC_TIME=French_Switzerland.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] vioplot_0.3.5 zoo_1.8-7     sm_2.2-5.6    knitr_1.30   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5       compiler_3.5.3   pillar_1.4.6     plyr_1.8.6      
##  [5] tools_3.5.3      digest_0.6.25    jsonlite_1.7.1   evaluate_0.14   
##  [9] lifecycle_0.2.0  tibble_3.0.1     gtable_0.3.0     lattice_0.20-38 
## [13] pkgconfig_2.0.3  rlang_0.4.7      yaml_2.2.1       xfun_0.18       
## [17] stringr_1.4.0    dplyr_1.0.2      generics_0.0.2   vctrs_0.3.4     
## [21] grid_3.5.3       tidyselect_1.1.0 qvalue_2.12.0    glue_1.4.2      
## [25] R6_2.4.1         tcltk_3.5.3      rmarkdown_2.4    farver_2.0.3    
## [29] ggplot2_3.3.2    purrr_0.3.4      reshape2_1.4.4   magrittr_1.5    
## [33] scales_1.1.1     ellipsis_0.3.1   htmltools_0.5.0  splines_3.5.3   
## [37] colorspace_1.4-1 labeling_0.3     stringi_1.4.6    munsell_0.5.0   
## [41] crayon_1.3.4